Basics

?Arithmetic
?Special

Linear Algebra

Vectors

The inner product

The inner product is also called the dot product or the scalar product. It is the sum of the item-by-item product.

> u <- rep(3,3)


> v <- 1:3


> u%*%v # the inner product
     [,1]
[1,]   18

The outer product

The outer product is also called the cross product or the vector product. It is a matrix resulting from the product of the elements of the two vectors.

> v <- rep(3,3)


> u <- 1:3


> u%o%v # The outer product
     [,1] [,2] [,3]
[1,]    3    3    3
[2,]    6    6    6
[3,]    9    9    9

Matrix Algebra

If you want to create a new matrix, one way is to use the matrix() function. You have to enter a vector of data, the number of rows and/or columns and finally you can specify if you want R to read your vector by row or by column (the default option) with byrow. You can also combine vectors using cbind() or rbind(). The dimension of a matrix can be obtained using the dim() function or alternatively nrow() and ncol().

> matrix(data = NA, nrow = 5, ncol = 5, byrow = T)


> matrix(data = 1:15, nrow = 5, ncol = 5, byrow = T)


> v1 <- 1:5


> v2 <- 5:1


> cbind(v1,v2)


> rbind(v1,v2)


> dim(X)


> nrow(X)


> ncol(X)

Some special matrix

The identity matrix has ones on the diagonal and zeros outside the diagonal.

• eye() (matlab)
• diag(1,nrow=10,ncol=10)
• diag(rep(1,10))

J matrix is full of ones

• ones() (matlab)

A matrix full of zeros

• zeros() (matlab)

> library(matlab)


> eye(3)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1


> ones(3)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1


> zeros(3) 
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

Diagonal matrix

> diag(3)

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

Upper triangular

> round(upper.tri(matrix(1, n, n))) 

for n=3
     [,1] [,2] [,3]
[1,]    0    1    1
[2,]    0    0    1
[3,]    0    0    0

If you also need the diagonal of one's 



> round(upper.tri(matrix(1, 3, 3), diag = TRUE))

      [,1] [,2] [,3]
[1,]    1    1    1
[2,]    0    1    1
[3,]    0    0    1

Lower triangular

Same as upper triangular but using lower.tri instead

• create an Hilbert matrix using hilbert() (fUtilities).

Matrix calculations

• compute a matrix multiplication X%*%Y.

> b <- matrix(nrow = 2, ncol = 2, c(1, 2, 3, 4))


> a <- matrix(nrow = 2, ncol = 2, c(1, 0, 0, -1))


> a
     [,1] [,2]
[1,]    1    0
[2,]    0   -1


> b
     [,1] [,2]
[1,]    1    3
[2,]    2    4


> a%*%b
     [,1] [,2]
[1,]    1    3
[2,]   -2   -4


> b%*%a
     [,1] [,2]
[1,]    1   -3
[2,]    2   -4

• compute the Kronecker product using %x% or kron() (fUtilities).

> M <- matrix(rep(2,4),nrow = 2) 


> M
     [,1] [,2]
[1,]    2    2
[2,]    2    2


> I <- eye(2) 


> I
     [,1] [,2]
[1,]    1    0
[2,]    0    1


> I %x% M 
     [,1] [,2] [,3] [,4]
[1,]    2    2    0    0
[2,]    2    2    0    0
[3,]    0    0    2    2
[4,]    0    0    2    2


> library(fUtilities)


> kron(I,M)
     [,1] [,2] [,3] [,4]
[1,]    2    2    0    0
[2,]    2    2    0    0
[3,]    0    0    2    2
[4,]    0    0    2    2

Matrix transposition

• Transpose the matrix

> t(M)
     [,1] [,2] [,3]
[1,]    1    0    1
[2,]    0    1    2
[3,]    0    0    1

The trace and determinant of a matrix

• compute the trace of a matrix using tr() (fUtilities)
• returns the rank of a matrix using rk() (fBasics:)

Matrix inversion

• Invert a matrix using solve() or inv() (fUtilities). We can also compute the generalized inverse using ginv() in the MASS package.

> M <- cbind(c(1,0,1),c(0,1,2),c(0,0,1))


> solve(M)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]   -1   -2    1


> solve(M)%*%M
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

Solving a linear equation

> m=matrix(nrow=2,ncol=2,c(1,-.8,1,.2))


> m
     [,1] [,2]
[1,]  1.0  1.0
[2,] -0.8  0.2


> 


> l=matrix(c(1.0+25.0/18,25.0/18.0))


> l
         [,1]
[1,] 2.388889
[2,] 1.388889


> 


> k=solve(m,l)


> k
           [,1]
[1,] -0.9111111
[2,]  3.3000000


> 


> m%*%k          #checking the answer
         [,1]
[1,] 2.388889
[2,] 1.388889
>

Eigenvalue, eigenvector and eigenspace

• Eigenvalues and eigenvectors

> eigen(M)
$values
[1] 1 1 1

$vectors
     [,1]          [,2]          [,3]
[1,]    0  2.220446e-16  0.000000e+00
[2,]    0  0.000000e+00  1.110223e-16
[3,]    1 -1.000000e+00 -1.000000e+00

Misc

• compute the norm of a matrix using norm() (fUtilities).
• check if a matrix is positive definite isPositiveDefinite() (fUtilities).
• make a matrix positive definite makePositiveDefinite() (fUtilities).

• computes row statistics and column statistics (fUtilities).
• extract the upper and the lower part of a matrix triang() and Triang() (fUtilities).
• See also the matrix, matlab, matrixcalc, matrixStats packages.

Analysis

Logarithm and Exponents

We have the power function 10^3 or 10**3 , the logarithm and the exponential log(2.71), log10(10),exp(1).

> 10^3 # exponent
[1] 1000


> 10**3 # exponent
[1] 1000


> exp(1) # exponential
[1] 2.718282


> log(2.71) # natural logarithm
[1] 0.9969486


> log10(1000) # base 10 logarithm
[1] 3


> log(1000,base = 10) # base 10 logarithm
[1] 3

Polynomial equations

To solve ${\displaystyle ax^{k}+bx^{k-1}+\cdots +n=0}$, where ${\displaystyle a,b,\dots ,n}$ are given numbers, use the command

> polyroot(c(n,...,b,a))

So, for example, to calculate the roots of the equation ${\displaystyle 2x^{2}-5x-3=0}$ one would do as follows:

> polyroot(c(-3,-5,2))
 [1] -0.5+0i  3.0-0i

and the solution can be read to be ${\displaystyle x_{1}=-0.5\land x_{2}=3}$.

See also polynom and multipol packages

Derivatives

Symbolic calculations

R can give the derivative of an expression. You need to convert your function as an expression using the expression() function. Otherwise you get an error message.

Here are some examples :

> D(expression(x^n),"x")
x^(n - 1) * n


> D(expression(exp(a*x)),"x")
exp(a * x) * a


> D(expression(1/x),"x")
-(1/x^2)


> D(expression(x^3),"x")
3 * x^2


> D(expression(pnorm(x)),"x")
dnorm(x)


> D(expression(dnorm(x)),"x")
-(x * dnorm(x))

Numerical approximation

• numDeriv package

Integration

R can perform one dimensional integration. For example we can integrate over the density of the normal distribution between ${\displaystyle -\infty }$ and ${\displaystyle +\infty }$

> integrate(dnorm,-Inf,Inf)
1 with absolute error < 9.4e-05


> integrate(dnorm,-1.96,1.96)
0.9500042 with absolute error < 1.0e-11


> integrate(dnorm,-1.64,1.64)
0.8989948 with absolute error < 6.8e-14
# we can also store the result in an object


> ci90 <- integrate(dnorm,-1.64,1.64)


> ci90$value
[1] 0.8989948


> integrate(dnorm,-1.64,1.64)$value
[1] 0.8989948

see the adapt package for multivariate integration.

> library(adapt)


> ?adapt


> ir2pi <- 1/sqrt(2*pi)


> fred <- function(z) { ir2pi^length(z) * exp(-0.5 * sum(z * z))}


> 


> adapt(2, lo = c(-5,-5), up = c(5,5), functn = fred)
       value       relerr       minpts       lenwrk        ifail 
    1.039222 0.0007911264          231           73            0 


> adapt(2, lo = c(-5,-5), up = c(5,5), functn = fred, eps = 1e-4)
       value       relerr       minpts       lenwrk        ifail 
    1.000237 1.653498e-05          655          143            0 


> adapt(2, lo = c(-5,-5), up = c(5,5), functn = fred, eps = 1e-6)
      value      relerr      minpts      lenwrk       ifail 
   1.000039 3.22439e-07        1719         283           0

• See also integrate.gh() in the ecoreg package.

Probability

• The number of combination of length k within n numbers : ${\displaystyle C_{n}^{k}}$

> choose(100, 5)
[1] 75287520

• Union and intersection

> union(1:10, 5:7)
[1]  1  2  3  4  5  6  7  8  9 10


> intersect(1:10, 5:7)
[1] 5 6 7

Arithmetics

The factorial function

factorial returns the factorial of an integer. This can also be computed using the prod() (product) applied to the vector of integers between 1 and the number of interest.

> factorial(3)
[1] 6


> prod(1:3)
[1] 6

Note that by convention ${\displaystyle 0!=1}$. factorial() returns 1 in 0. This is not the case with the prod() functions.

> factorial(0)
[1] 1


> prod(0)
[1] 0

Factorial numbers can be very large and cannot be computed for high values.

> factorial(170)
[1] 7.257416e+306


> factorial(171)
[1] Inf
Message d'avis :
In factorial(171) : value out of range in 'gammafn'

The modulo function and euclidian division

• Modulo and integer division (i.e. euclidean division)

> 5%%2
[1] 1
>5%/%2
[1] 2

Note: R is affected by the problem with non integer numbers and euclidian divisions.

> .5%/%.1 # we get 4 instead of 5
[1] 4


> .5%%.1 # we get .1 instead of 0
[1] 0.1

Geometry

• pi the constant
• cos(), sin(), tan() the trigonometric functions.

Symbolic calculus

rSymPy (rsympy) provides sympy (link) functions in R.

If you want to do more symbolic calculus, see Maxima^[1], SAGE^[2], Mathematica^[3]

See also

The following command gives help on special mathematical functions related to the beta and gamma functions.

?Special

References